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Abstract 

We report the results of extensive numerical experiments on the Kuramoto-Sivashinsky equation 
in the strongly chaotic regime as the viscosity parameter is decreased and increasingly more linearly 
unstable modes enter the dynamics. General initial conditions are used and evolving states do not 
assume odd-parity. A large number of numerical experiments are employed in order to obtain 
quantitative characteristics of the dynamics. We report on different routes to chaos and provide 
numerical evidence and construction of strange attractors with self-similar characteristics. As the 
“viscosity” parameter decreases the dynamics becomes increasingly more complicated and chaotic. 
In particular it is found that regular behavior in the form of steady state or steady state traveling 
waves is supported amidst the time-dependent and irregular motions. We show that multimodal 
steady states emerge and are supported on decreasing windows in parameter space. In addition we 
invoke a self-similarity property of the equation, to show that these profiles are obtainable from 
global fixed point attractors of the Kuramoto-Sivashinsky equation at much larger values of the 
viscosity. 
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1 Introduction 


The Kuramoto-Sivashinsky equation is one of the simplest one-dimensional PDE’s which exhibits 
complex dynamical behavior. As an evolution equation it arises in a number of applications includ- 
ing concentration waves and plasma physics ([4], [20], [21], [22]), flame propagation and reaction 
diffusion combustion dynamics ([30], [31]), free surface film-flows ([2], [14], [29]) and two-phase 
flows in cylindrical or plane geometries [7], [11], [25]. [35], [36]. The equation can be written as 


L } t + UV X + U XX + UxTXX — 0 . 


r 1 L ^ 

x € [ 2 1 2 per ' 


( 1 ) 


noting that the solution remains periodic of period L if the initial data are periodic. 

Equation (1) can be normalized to 27r-periodic domains by the transformations X\ = xx/L, 
U — xu/L and t\ = x 2 t) L 2 , which on dropping the subscripts 1 yields the scaled KS equation 


Ut + uu x 4- ^ xx 4" V^xxxx — 0: 

(x,t) € 1R 1 x ]R + , (2) 

u(x , t ) = u(x 4- 2?r, J), 

where v = £2 > 0 represents a dimensionless wavelength of u. Physically, then, the limit v — ► 0 
corresponds to the waves in (2) becoming infinitely long. It is easily established by linearization of 
(2) that for a given v the first mod(i/“ 2 ) Fourier modes are linearly unstable and grow exponen- 
tially. It is well-known, from numerical experiments for instance (see below), that solutions to (2) 
become increasingly irregular as v is decreased. 

Without loss of generality the study of solutions with zero spatial mean is sufficient due to 
Galilean invariance. An analytical study of the KS equation was carried out by Nicolaenko et al. 
[24] (referred to as [NST]), where it is shown that if the initial data are in L 2 and are of odd-paritv 
(i.e. antisymmetric about x — 0 for (1) or x = x for (2)), then the solutions remain in L 2 for all 
time and there is a globally attracting set also bounded in L 2 . The bound depends on the value of 
L or v and the estimate of [24] gives the following upper bound: 

lim sup||u(*, t)| I 2 < const • u~ 3 ^ 2 . (3) 


The corresponding bound for (1) is 0(L 5 / 2 ). Such bounds of the L 2 norm are useful in determining 
estimates of the Hausdorff dimension, d//, of the universal attractor. In fact working with equation 
(1), it has been shown by Temam [34], that if the bound for the L 2 norm is 0{L Q ) then dfj is 
proportional to Z,( 7 4- 2of )/ 8 and so the [NST] result yields the estimate L 13 / 8 for d//; we note that 
the conjectured best bound is O(L) but a proof is not available yet. 

The analysis in [NST] was based on the odd-parity of the initial data. This restriction has 
been removed recently by several authors. Goodman [13] considers general smooth initial data and 
obtains a bound of the same form as (3). This bound is improved further, however, by ITy&shenko 
[17] and independently using a more classical approach by Collet et al. [5], and is 

lim sup||u(*, t)| I 2 < const • j/“ 13 / 10 . (4) 

t — ►oo 


Again the analogous bound for (1) is 0(L s ^ b ) improving the estimate of the Hausdorff dimension 
to X 51 / 40 bringing it closer to the conjectured best bound. We note that the analysis in [5] and [17] 
applies also to odd-parity initial conditions and so is an improvement of the method in [NST]. 

Analyticity of the solution has also been proved by Collet et al. [6]. In fact, considering equation 
(1), the result proved is that for large X the function U(x,t) is analytic in a: in a strip of width 
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(3 > const • L 16//25 around the real axis, which in turn implies that the high frequency part of the 
spectrum has the form 

\U n (t)\ « 0(exp(-const • L~ ie f 25 q\n\), (5) 

where U n (t) is the nth Fourier coefficient of U(x,t) and q = 27 r/X. By a series of numerical 
experiments a much stronger result that (5) has been conjectured in [6], namely that there exists 
a j3 > 0 independent of L such that the solutions of (1) satisfy 

lim sup ^ | e /?<? l n l{7n(0| < const., (6) 

neZ 

and the numerical work gives (3 ^ 3.50. These results are extremely useful in guiding our numerical 
work as is explained in Section 2. 

There are several computational and analytical studies dealing with the construction of approx- 
imate inertial manifolds of KS. Kirby [19] develops a Galerkin approximation based on Sobolev 
eigenfunctions; Chen [3] constructs approximate inertial manifolds by broadening the gaps in the 
spectrum in order to obtain low dimensional behavior; Foias et al. [10] construct fully discrete 
nonlinear Galerkin schemes based on approximate inertial manifolds of KS; Robinson [28] uses the 
bounds of [5] and [17] to produce a new estimate of the dimension of the globally absorbing set of 
KS. Additional references may be found in these articles. 

A significant amount of numerical work preceded the theoretical results outlined above. The 
numerical experiments revealed a wealth of interesting nonlinear phenomena and in particular 
highly complex dynamics including chaotic trajectories as the dissipation parameter v decreases. 
Earlier works include the computations of Cohen et ah [4], Sivashinsky and Michelson [32], Aimar 
[1], and Manneville [23]. Systematic explorations of phase space were carried out by Hyman and 
Nicolaenko [15] and Hyman et al. [16], Papageorgiou and Smyrlis [27], [33] and Coward et al. [7] 
who report on many features of the dynamics, Kevrekidis et al. [18] computed the bifurcation 
diagram for relatively large values of */, using a bifurcation algorithm. For smaller values of v 
unsteady phenomena set in; it was first found by Papageorgiou and Smyrlis [27] and also Smyrlis 
and Papageorgiou [33], that in the case of odd-parity initial conditions there is a period-doubling 
route to chaos. Extensive numerical solutions were employed to provide strong evidence that the 
route to chaos is according to the Feigenbaum scenario ([8], [9]); in addition the universal constants 
computed by Feigenbaum for one-dimensional nonlinear non-invertible maps, were computed from 
our numerical data with two digit accuracy. In the sequel we term chaotic dynamics just beyond 
the accumulation point, in parameter space, as Feigenbaum chaos. 

In this article we present the results of extensive numerical computations for general initial 
conditions. Comparisons with previous studies (e.g. [15], [16], [18]) have been made where there is 
an overlap with full agreement. Our particular interest is in lower values of v where the dynamics 
gets increasingly complicated. 

2 Numerical methods 

2.1 Approximation of the solution 

The assumption that the initial data is spatially periodic of period 2x, allows us to represent the 
solution u(x,t) of the Kuramoto- Sivashinsky equation (from now on KS) as : 

u(x,t)= cos kx + (3k(t) sin kx) + ao(t) (7) 

k> i 
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The conservative nature of KS implies that a >o(t) remains constant, i.e. 


doo(t) 

dt 


1 d ' 2 * 
2 tt dt 

= — — / 2 '( 
27t Jo 


y" 2 7T ^ /*2?r 

/ = — / Ut(x,f)<i 

Jo 27T Jo 


UU X — Wjj — f— r^U xxxx )ofx — 0. 


One could easily observe that whenever u(x,C is a solution then so is u(x - cEf) + c, which allows 
us to assume ao(t) = 0 for simplicity. Replacing u(x,t) by its Fourier representation in the PDE 
we obtain 


V>t 4~ UU X 4" ^xx 4” — 


where 


= - (A: 2 - vk 4 )a k - /U) cos /ex + (j3' k - ( k 2 — vk 4 )j3 k — B k ) sin kx) 

k> 1 


= - 


* 

2 


^ , Q-mPn 4" g 

m+n=/: 


Ot n j3 m ) 

77{ — 71— k 


and 

k k 

Bk — ~ ^ ^ ( ft m Q n ~~ fim&n) 4“ “ } . ( a m Q n 4“ Pmfin) 

m + n—k m—n=k 

Thus finally KS equation becomes equivalent to an infinite dimensional system of ODEs : 

OLk ~ + Ak 

flk = ^kfik 4" B k 

where k 6 IN and = k 2 — uk 4 w r ith 

Afc = ^(01,02, --'iPufo, ••■) an d Bk — 5jt(ai, «2? ^2, ...) 

In the case of the Generated AJ/ramofo-Siuas/nns/ey equation 


( 8 ) 


(9) 


(10) 

( 11 ) 


+ uu x + u xx + ^u xxxx 4- Pu = 0 


( 12 ) 


where P is a dispersive kernel, see Papageorgiou et al [25], the resulting infinite dimensional system 
of ODEs is quite similar 

o]t = + <5/3* 4- A k (13) 

P k = ^ikAfe “ <5a* 4" i?* (14) 


where £ is defined by 


( Vu)(k ) = t£u(fc) 


(15) 


It is already established by Hyashenko [17], Collet et al [5], Goodman [13], that all the spatial 
Sobolev norms of the solution of KS equation for arbitrary initial conditions, remain bounded for 
all times, which implies that 

«* = hm sup|a* 2 + 0 k 2 \ 1/2 

f — + OC 
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decays faster than any algebraic rate. Collet et al. [6] have already proved that the decay is 
exponential. Furthermore they present numerical evidence to support a conjecture that the number 
of determining Fourier coefficients is proportional to r/” 1 / 2 . 

Such numerical evidence justifies approximation of the solution of the KS by truncation of 
higher frequencies. The size of the truncation can be determined by the Collet et al. [6] numerical 
study of the decay rate of the Fourier coefficients, but a first estimate of this is obtainable from 
the argument that follows. The KS equation in Fourier space can be written as : 

■^ru(k, t) — \ku(k, t) — ( uu x )(k y t) 
at 

or 

^w(M) = A*u(M) - *fc(^« 2 )(M) 

We also have : 

Fact : If fji is positive and f(t) is in Z°°[0, -hoc) then the solution of x* + fix = f(t) satisfies: 

lim sup |x(t)| < 1 1 1/| |oo (18) 

*— ►+ oo // 


(16) 

(17) 


Using the above result in combination with (17), we observe that, if we had, even a rather rough 
estimate of the size of u and its norms, we could obtain an idea of the number of Fourier coefficients 
with significant numerical contribution to the approximate solution of the KS equation. 

Having used various numerical schemes, we obtained sufficiently good estimates of the size of 
u, as well as of u 2 , uu x and of \u\tf. Thus we subsequently obtained an estimate of the number of 
numerically significant Fourier coefficients. It should be mentioned here that the existing analytical 
results alone can not be of much practical value in determining the truncation size. 

Currently we use a finite difference scheme of variable time step to integrate the truncated 
system of ODEs 

Oik = ^kOtk + Ak(ai,<X2i"‘'><XN,0U02,*"'>0N) ( 19 ) 

Pk — ^kfik + &k ( Q '1^ Q 2> (20) 

where k = 1,2 , N and 

k k ~ l k N ~ k 

A k = — 2 ^ "j Omftk—m T X y v ( G m0m+k ~ Omflm+k ) (21) 

Z m= 1 1 m= 1 


k k N ~ k 

^ ^ (o m Ok-m ~ 0m0k—m) + ^ ^ ^ T Pmfim+k ) (22 ) 

m— 1 m= 1 


Our scheme turns out to be extremely stable due to the form of the linear part of the right hand 
side of the system. We emphasize that even though A* is bounded from above, with a mild upper 
bound namely 


A* < 


4i/ ’ 
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for all k £ IN, it could have a negative lower bound very large in absolute in our truncated system. 
For example, if N = 32, v - .01 then « -9.5 x 10 3 . Thus if At is such that 

\l N \At = 0(1) 

then the relative local numerical error for a^v, (3n is also 0(1), since at each time step the linear 
semigroup multiplies a/v, Pn by exp(lj^At). We should also note here that in the nonlinear part of 
the system of ODEs corresponding to high frequencies in our truncated system, i.e., B^\ most 

of the contribution comes from the low frequency coefficients. This implies that higher frequency 
coefficients are slaved by the low frequency ones due to the scheme; this is a distinguishing feature 
of dissipative infinite dimensional dynamical systems and is thus inevitable. Still we wish to allow 
the high frequencies as much freedom as possible in order to enable potential individual behavior. 
In achieving this we incorporate a modification suitable for stiff systems of ODEs. Stiffness is 
caused here by the big variation of the values of A^ and thus the time step in our modified scheme 
varies according to the size of A^. Nevertheless the relative local error in every individual term also 
varies, in our scheme, according to its contribution to the solution. 

Altogether the local numerical error, due to the approximation of the solution of the system of 
ODEs, is reduced to be insignificant compared to the one caused by the truncation alone. It should 
be clarified at this point that even with coarser time step subdivisions than the ones in current use, 
the relative local error to the solution profile due to the high frequencies is satisfactory. High order 
accuracy, i.e. the local relative error being kept below the level of 10~ 6 . is employed only for the 
sake of studying high order frequencies and their long time influence on the nature of the global 
attractor. 


2.2 Numerical Experiments 

We have carried out a very large number of numerical experiments, testing more than 800 values 
of the dissipation parameter v ranging from* v = .99999999 to v = .002. In representative values 
of our dissipative parameter v we carried out several runs to make sure that the behavior of the 
attractor corresponding to that v does not change with more Fourier coefficients of finer time 
step. This is necessary since our major concern is not just the local error but faithful reproduction 
of the long time behavior. The size of the truncation ranged from N = 4 to N = 512 making 
computation in the latter case rather slow. Most important of all is that many cases had to be 
followed for very long time in order to achieve convergence to the corresponding attractor and/or 
to accurately classify the attractor. For example in the case v — .1212267996068, the attractor is 

*If v > 1 then 

1 [ 2 * 

lim u(x,t) — wo = — / uo(x)dx, (23) 

2t r J Q 

since KS equation implies that 

2~dt i u (’ , *)l^ = DU — *d u **(’> *)h 

and for u £ I 2 [0,27t], 27r-periodic, with f** u(x)dx — 0 we have the following inequalities: 

|«h < \u x \ 2 < \u xx \ 2 

On the other hand if v < 0 the problem is ill-posed. 
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periodic with 2 14 distinct maxima and as many minima in the L 2 norm as a function of time. This 
particular value of v follows previous ones corresponding to 13 period doubling bifurcations. The 
period is estimated to be T « 25573, which corresponds to approximately 5.1 x 10' time steps and 
each step corresponding to substeps. Satisfactory convergence to the periodic attractor occurred 
after 10 9 time steps. 

In representative cases, of different values of v, we have used our method to estimate the error 
caused by the truncation. The relative error, with respect to the L 2 and H 1 norms ranged from 
10 -6 to numerically insignificant. Though 10 -6 might look big, we should point out that our 
study’s major target was not to run the most accurate algorithm - which we did in selected cases. 
Instead it was to carry out many and long runs, and keep sufficient accuracy, in order to reproduce 
the nature of the corresponding attractor. Classification of the nature of the attractors at each 
different value of v made it possible to determine the windows in the parameter v space. Most of 
the runs had to be carried out in order to accurately determine the limits of the windows of v. and 
furthermore determine whether the transition between the two kinds of attractors was smooth or 
not. A case of a smooth transition is a bifurcation; A period-doubling is one such as well as a case of 
eigenvalues bifurcations of the Jacobian of the flux of our dynamical system which cause transition 
from stationary /traveling attractors to periodic or chaotic ones. A non-smooth transition is a case 
when we have for a i/-interval, [j/j, 1 / 2 } say, coexistence of two or more invariant locally attracting 
sets. While one of them was attracting most of the initial data, as v varies, another one becomes 
suddenly more likely to attract most initial data. For example, in the case of v = .232, we have 
coexistence of a unimodal traveling attractor and of one with periodically appearing homoclinic 
bursts. The second one is more attracting. Even though the exact bifurcation values of u depend 
on the accuracy of the numerical scheme, the nature of the bifurcation does not depend on the 
accuracy provided that the method is of sufficient accuracy since the system has finite degrees of 
freedom. This observation forces us to use higher accuracy when we seek better approximation of 
the bifurcation values of u. 


2.3 Characterization of the attractors 


As v varies the long time behavior of the solution does also, in some cases smoothly and in others 
not. Different values of v correspond to stationary (unimodal, bimodal, trimodal,...), traveling, 
periodic and chaotic attractors. Stationary attractors are easy to observe, and their stability can 
be studied with standard methods. It suffices to check the eigenvalues of the Jacobian of the 
truncated system at the stationary profile. Similarly in the case of traveling profiles, where one 
could easily determine its speed of propagation : 


c = X7 Are < ’ „? + ,)? > + 0( A * ' 


(24) 


where 


(«!,&) = (Qi(t- At),fh{t- At)) 


(a 1 J i ) = (a 1 (t + At),p 1 (t + At)) 


In the case of periodic attractors more sophisticated numerical methods are required for their 
quantification. First suspicion of a periodic attractor comes from the energy versus time plot - 
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E(t) = \u(t, -)\ L 2 . The evidence becomes more convincing once we obtain the phase plane of the 
energy plots - i.e. E(t) versus dE(t)/dt, where the values of dE(t)/dt are accurately obtainable 
using a suitable interpolation of several points (4? E{ik)) . Periodic attractors correspond to closed 
phase curves and the number of minima/maxima is the number of the points where the phase 
curves intersect the E — axis , i.e. dE(t)/dt, = 0 with d 2 E(t)/dt 2 > 0 or < 0 respectively. 

Nevertheless we wish to have more accurate quantitative data such as the exact period as well as 
the exact values of the local extrema, in order to accurately determine period doubling bifurcations 
and study the potential fractal nature of the return map§ of the extrema. In order to achieve that, 
we have used a method of calculation of the extrema by an optimized polynomial interpolation over 
a sufficiently large number of points (t*, £(<*)) . The polynomial of the interpolation is obtained 
by a suitably weighted least squares method over consecutive points, many more than the degree 
of the polynomial sought for. Weights in the least squares approximation, reduce the effects of 
higher order polynomial terms and the round-off error caused bv the computation of the quantity 
used in the interpolation. It is of critical importance that all points of the interpolation lie in the 
convex or concave region, when we look for minima or maxima respectively. This method, in the 
case of test functions, gave us accuracy higher than the machine precision, nevertheless, as it could 
be shown analytically, the precision in the estimation of local extrema, depend on the time step 
and the degree of the interpolation polynomial. Weights are introduced to make the points which 
are closer to the local extremum have a higher contribution to the interpolation. 

The original motive in developing the interpolation algorithm was to accurately determine the 
period doubling bifurcations and establish the appearance of the Feigenbaum universal constants. 
(See Smyrlis and Papageorgiou [33].) This method turned out to be extremely helpful in deter- 
mining features such as the fractal nature of attractors, by detecting such behavior in the return 
map of the energy minima, for example. Most successfully it has enabled us to detect and classify 
the quasiperiodicity of attractors. We should mention here that the same conclusions, concerning 
periodic attractors are also obtainable if instead of the L 2 norm and its extrema we had used for 
their detection, other quantities such as the H 1 norm and the Fourier coefficients. However there 
is an exception here, the case of traveling/periodic attractors, where instead of 

u(x, t 4- T) = u(x,t) for all {x,t) G IR x IR + 


we have 


u(x,t + T) = u(x — cT.t) for some c € IR 

Thus E(t) has period T, but ii(k.t) is in general a quasiperiodic function unless cT/27r is a rational 
number. Fourier coefficients being quasiperiodic produce very beautiful pictures - ak(t) versus 
3k(t )- The speed of propagation and period could also be accurately estimated using a suitable 
interpolation. 

In the classification of chaotic attractors and transitions to chaos scenarios, then, the extrema 
of E(t) and their return maps played the key role. We should mention here that if f(t) is a 
quasiperiodic function, for example 

f(t ) = cos t + cos(\/2 1), 

then, although there is not much to observe in the graph, one could observe a lot in the return 
map. If (mjfc)fce ]\ is the sequence of the local minima, in the order they appear for t > 0, then their 

s Let a„, n £ IN be a sequence of real numbers, then the set of points (a n ,an+i), n € IN, which is a subset of IR'. 
is called the Return Map of (a„)neiN- In the case of a function / : [0,oo] — IR the return map of / is the return map 
of the sequence of the local minima (or maxima or extrema) of /. 



return map, lies on a dosed curve. In general if /i(i), f 2 (t) are periodic functions with irrational 
frequency ratios, then the return map of f(t) = f\{t) + f 2 (t) lies on one or more lines. This is 
generalized to three or more irrational frequencies. 

Visualization of the return map of the local minima (or maxima) of the energy (or other quanti- 
ties such us the H 1 norm which magnifies the contribution of high frequencies) enables us to classify 
the attractors, not only in the quasiperiodic case, but also in more complicated chaotic ones. 

In certain cases the return map plot has a fractal nature. Foliations are observable when one 
plots successive magnifications of the graphs. We are in the process of estimating numerically the 
HausdorfF dimension of the graph corresponding to the return map in order to establish the rate of 
growth of the HausdorfF dimension as v decreases, particularly in those chaotic cases where there is 
no recognizable pattern. In such studies it is imperative that the local extrema are computed with 
a high accuracy, in particular if magnifications in the return map of several orders of magnitude 
are to be carried out. 

3 Numerical results 

3.1 Dynamics for “large” v\ 0.09 < v < 1 

For values of v larger than about 0.1 a fairly complete picture of the dynamics has been given 
by numerous authors ([15], [16], [18]). Steady state or traveling waves, including their bifurcation 
diagrams are computable by well- developed bifurcation algorithms, see for example ([18]). Our 
interest is on the spatio-temporal evolution in chaotic attractors which led to the development of 
the accurate code tracking the time evolution of solutions, described in Section 2. For completeness, 
however, and to set the stage for the results that follow we give a brief description of the dynamics 
corresponding to “large” v . We note that as v decreases computed attractors are not necessarily 
global ones; our numerical results are based on large time solutions of general initial value problems 
and the numerical solutions given here belong to the most strongly attracting ones in the sense 
that they are reproduced by fairly different initial conditions (e.g. a sinusoidal initial condition or 
one having a white noise spectrum). 

A schematic of the dominant attractors in this range is given in Figure 1. The interval v > 1 
is not included in Figure 1 since it is trivial in the sense that a uniform steady state is achieved 
here which is equal to the spatial mean of the initial data. The figure is drawn to scale in order 
to give the relative sizes of attractors. The letters A-F are used to identify the different inter- 
vals which are briefly described below. According to extensive computations, there is an interval 
(.12116, .1212267996064) which cannot be distinguished on Figure 1 due to its small size - w r e call 
this interval E\ and describe it below. 

Interval A: 0.307602 < v < 1. 

For 0.30765 < v < 1 solutions get attracted to a 27r-periodic steady state (unimodal fixed point ). 
The energy, or the spatial X 2 -norm of the solution, defined as E(t ) = u 2 (x y t)dx increases as a 
function of v in this interval and attains a maximum at v w 0.59094. 

Interval B: 0.232491 < v < 0.307601. 

For 0.232491 < v < .307601 a 27r-periodic (unimodal) steady-state traveling wave is found - the 
wavespeeds increase monotonically from zero near the upper edge of the window to approximately 
0.956 at v = 0.232491. 
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Interval C: 0.17735 < v < 0.23249. 

A bifurcation occurs just below v = 0.232491, and in the range 0.17735 < v < 0.23249 periodic 
homoclinic bursts (PHB) emerge as the most attracting states. The PUB solutions found in this 
window have Z 2 -norms which are constant except for periodically occurring homoclinic bursts. The 
time intervals between bursts vary from approximately 21 — 30 time units for v between 0.23249 
and 0.2, and then increase dramatically as v decreases further; for example the time between bursts 
for v = 0.179,0.178,0.1775,0.1774,0.17735 are approximately 200,450,1700,4500 and 18000 time 
units respectively. 

Representative results are given in Figures 2 and 3 for v — 0.2. Figure 2 shows the evolution 
of E(t) after the solution has settled into an attractor - a range between 750 and 1000 time units 
is shown. The energy is constant with homoclinic bursts occurring at approximately every 30 time 
units. The phase plane of E(t) provides strong evidence that the bursts are periodic. The solution 
spends most of its time at the rightmost point in the phase plane diagram and every 30 time units it 
loops around in less than 5 time units, with the process repeating periodically. The profile between 
bursts is steady for 30 time units and is in fact a bimodal fixed point, that is only even Fourier 
modes are non-zero. It has been confirmed that the bimodal profile during bursts is obtainable 
from the unimodal global fixed point attractor at u = 0.8 by use of the property (26) discussed 
in Section 2. We can surmise, therefore, that the bimodal fixed point attractor is unstable and 
overlaps with a time-periodic attractor. 

The evolution of E{t) does not provide details about individual Fourier coefficients since it is an 
integral quantity. The time history of the first two harmonics is given in Figure 3. The coefficients 
of cos(x), sin(x) and cos(2x), sin(2x) are plotted on the same scale so that their evolution can be 
followed concurrently. It is evident from the primary Fourier mode plots that the amplitudes are 
zero except during bursts which appear as sharp peaks. The peaks alternate from small size to larger 
size during each burst; there also seems to be an energy exchange between the two amplitudes in 
that during each burst the cos(x) and sin(x) amplitudes attain alternatively large and small peaks. 
The sign of individual amplitudes during bursts does not appear to have a recognizable pattern 
but seems to be random. 

The situation is much more regular for the 2nd Fourier mode. The time evolution of the 
amplitudes of cos(2x) and sin(2x) are now in phase with bursts appearing as sharp jumps connecting 
equal but opposite steady states. Similar results persist for higher modes, the main difference being 
a marked decrease in amplitude. 

Interval D: 0.131815 < v_ < 0.1773. 

The time increase between bursts heralds the appearance of a new fixed point attractor. It is 
found that in the range 0.131815 < v < 0.1773 solutions get attracted to bimodal steady states, 
i.e. 7 r-periodic stationary profiles. According to the property (26) these solutions are obtainable 
from the unimodal solutions in the subinterval [0.5272,0.7092] of A by choosing p = 2. 

Intervals E, 0.1212267996068 < v < 0.13181, and E 2 : 0.12116 < v_ < 0.1212267996064. 

The bimodal steady states loose stability and a time periodic attractor is entered. The energy 
E{t) has an approximate period of 0.93862 time units at the beginning of the window, v = 0.13181. 
The period increases monotonically as u is decreased and the fitst period doubling takes place in 
the interval .12175 < v < .121752. This pattern is repeated, i.e. a monotonic increase in the 
period and then a period doubling, and the following eleven successive period doublings were found 
in each of the open intervals (.1213145, .121315), (.12124493, .12124495). (.12123065, .1212307), 
(.12122762, .12122775), (.121226975, .12122698), (.121226835, .12122684), (.121226805, .121226808), 
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(.1212268012, .1212268014), (.1212267998, .1212268), (.12122679968, .1212267997) and 
(.121226799622, .121226799625). The period by the twelveth period doubling is T a 12780 time 
units. Figure 4 shows the phase planes of the energy for representative cases capturing the first 
five period-doublings. The values of v are given on the Figure as well as the number of minima 
of the energy - there are 2 P minima where p is the number of period doublings that have taken 
place. It has been confirmed, but is not included here, that the cascade follows the Feigenbaum 
scenario ([8], [9] ) with //-intervals between successive period doublings decreasing geometrically by 
a universal factor 6 = 4.6692016.... The data analysis and methods used are identical to those 
of Papageorgiou and Smyrlis [27] and Smyrlis and Papageorgiou [33] who considered odd-paritv 
solutions of the KS equation. The second Feigenbaum universal constant which pertains to the 
geometry of the attractor is also accurately supported by those results. By analogy with the theory 
of non-invertible one-dimensional maps, coupled with the numerical evidence given here, it is rea- 
sonable to assume that an accumulation point exists, Vo say, below which the solution is chaotic; 
this of course is hard to prove since we are dealing with differential systems of large dimensions (24 
coupled equations are used in this window). 

After extensive computations we are able to identify an interval E\ = [0.12116, 0.1212267996064] 
which we describe as supporting Feigenbaum chaos. Evidence for this comes from the study of 
return maps of the energy minima (see Section 2 for a description of the data analysis tools used 
here) which give “strange attractors” with as many as 2 13 pieces at the beginning of interval E\. 
The number of pieces decrease and when v = 0.12118 the attractor consists of only one piece. A 
picture of a strange attractor constructed from the return map, is given in Figure 5(i)-(vii) which 
corresponds to v = 0.12122679907. The return map of a large number of energy minima beyond 
transient behavior is shown in Figure 5(i). There are 32,034 points in this plot and each point is 
represented by a dot. The points lie on a regular looking object and since each point gets mapped 
to another point in the plot by the action of the dynamical system, the pieces traced out mimic the 
analogous one-dimensional non-invertible map representation of the flow; the equivalent plot for 
the logistic map would trace out parts of a quadratic function. Interestingly, Figure 5(i) resembles 
the quadratic map but the differences have not been quantified. The similarity property of strange 
attractors is however exhibited in Figure 5. Starting from Figure 5(i), successive enlargements are 
made of the pieces of the attractor which are marked by the letter V. The scale of the enlargement 
can be followed by noting the horizontal scale of successive diagrams. The first enlargement, Figure 
5(ii), is of what appear as two pieces in 5(i) near the minimum of the plot. On enlargement a picture 
which is clearly similar to the complete attractor in 5(i) emerges. The two pieces near V in 5(ii) 
are enlarged to give 5(iii); again similarity with the previous plots is remarkable considering that an 
amplification of approximately a factor of 50 has taken place. The one piece near V in 5(iii) is now 
enlarged to give the three clusters in 5(iv). The rightmost marked piece is enlarged to give 5(v) and 
its rightmost piece is enlarged to give 5(vi). Finally an enlargement of the leftmost marked cluster 
of 5(vi) gives the final plot Figure 5(vii) which contains just four points. The successive plots 5(iv)- 
(vii) provide additional strong evidence of the self-similarity of the attractor. We note that such 
conclusions are only possible if a large and accurate data base is available - for instance in going 
from Figure 5(i) to Figure 5(vii), an amplification factor of approximately 3 x 10 6 has taken place 
along the horizontal axis. In the absence of a large and accurate data set of minima (see Section 2 
for the interpolation algorithms), the number of return map “iterations” which produce the strange 
attractor would be inadequate for sufficient levels of enlargement that exhibit self- similarity. In 
addition loss of accuracy in the minima estimation would produce a noisy data set which is of little 
value in exhibiting self- similarity. 
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Interval F: 0.09051 < v < 0.121155. 

Near the upper edge of this window, 0.121155 < u < 0.12105, solutions are attracted to ho- 
moclinic bursts, with the bursts being chaotic but with a low dimensional attractor. The bursts 
are spaced apart at roughly equal times and in between them the profile is bimodal and is obtain- 
able from the equivalent unimodal steady state in interval .4 according to the self-similar property 
(26). Typical results are given for v = 0.12115. The homoclinic bursting is shown in the evo- 
lution of the energy and its corresponding phase-plane, Figures 6; the phase-plane indicates the 
non-periodicity of the bursts. It has been confirmed that the profile in between bursts is given by 
u(x\ 0.12115) = 2U{2x: 0.4846) where U(y: v), 0 < y < 2n is the steady-state profile in interval A. 
Further information about the chaotic bursts is given in Figures 7(a)-(b). Figure 7(a) gives the 
evolution with time of the critical points of E{t), that is its local maxima and local minima. For 
a given time interval, there is a discrete set of energy maxima and minima and each of these is 
represented by a dot in the plot. The longer and relatively regularly spaced gaps with no points 
in the plot denote time intervals between bursts while darker segments correspond to the bursts. 
The irregularity of the positioning of the points indicates the chaotic aspect of the bursts and this 
can be seen more clearly in the return map of the energy minima in Figure 7(b) which suggests 
a strange attractor with foliations. The bimodal steady state in-between bursts is at the point 
(10.7415876, 10.7415876) on the return map diagram, which is marked with an open circle. Succes- 
sive enlargements of the attractor are given in Figure 8(i)-(vi). Again the point corresponding to 
the inter-bursting bimodal steady state is marked with an open circle and five successive enlarge- 
ments are made in its neighborhood. The similarity of the attractor is discernible, particularly if 
we note the horizontal extent of the attractor given in each plot in Figures (ii)-(vi); the horizontal 
extents are approximately 0.6. 0.15, 0.011, 0.001 and 0.000005 for Figures 8(ii)-(vi) respectively. 
An enlargement of the order 10 5 has taken place in going from Figure 8(ii) to Figure 8(vi) and 
the positioning of the dots, that is the geometry of the attractor, is remarkably similar; from a 
numerical point of view, we note that such details can only be captured if the energy minima are 
computed accurately to at least seven decimal places - this is another reason for the high order 
interpolation schemes used in constructing the minima from a set of energy values at discrete time 
intervals no smaller than the time-step of the integrator. 

As the value of v is decreased below approximately 0.121045 the liomoclinic bursts take place 
more regularly in time and they are almost periodic. A more careful scrutiny of the data from the 
representative case v = 0.1 reveals that the flow is chaotic, and the bursting phenomenon has a low 
dimensional attractor. Once again the unstable quiescent states between bursts give profiles w^hich 
are bimodal and given by u(x;0.1) = 2f7(2x; 0.4). A plot of E(t) together with a magnification of 
a burst is shown in Figure 9. The bursts are almost periodic and take place every 33 time units, 
approximately; the energy phase-plane has been shown not to be a closed curve of index 2 as in 
the homoclinic burst solutions found in interval C described above. Figure 10 shows the return 
map obtained from the energy minima. The point marked P corresponds to the steady state in 
between bursts and as the magnification shows, the attractor in the neighborhood of this fixed 
point consists mainly of two lines meeting at an angle smaller than x/2. Additional information 
about the bursting is given in Figure 11. The Fourier amplitudes of the first four modes are plotted 
against time, that is we plot Jaf + bf for i = 1,4 where u(x,t) = E^ =1 (a n cos(nx) + b n sm(nx)) 
is the Fourier-Galerkin representation of the solution. In between bursts, only the even modes are 
non-zero and as seen from the plots the amplitude of the fourth mode is smaller than that of the 
second one. During a burst, however, the even modes loose energy to the odd modes as can be 
seen clearly from the Figure; the amplitude of even modes always decreases (except right at the 
end of the burst) and at the same time the amplitudes of odd modes become non-zero for the 
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duration of the burst but vanish again after the homoclinic burst is finished. This behavior is very 
characteristic of the global energy transfer mechanism present in the KS equation. 

3.2 Dynamics for lower values of v: v < 0.09. 

As v decreases the dynamics become more complicated. A larger number of determining modes is 
required and co-existence of attractors is more likely. From a numerical viewpoint the computations 
become more expensive and accuracy criteria more stringent. 

A large number of numerical experiments have been carried out over several years and a sum- 
mary of the attractors is given in Figure 12. The Figure represents the i'-line (not drawn to scale) 
along with the most strongly attracting large time solutions found by us; note that the various 
attractors given in Figure 12 are not necessarily global ones but the window boundaries given 
delineate the regions where an attractor loses stability. The type of attractor was determined by 
applying the data analysis described earlier. 

In interval E described earlier, we provided evidence of a period-doubling route to chaos, and 
also a route to chaos through homoclinic bursting as in interval F above. A third route to chaos 
is via a quasi-periodic attractor and this route is supported by the first three windows in Figure 
12. A time periodic attractor looses stability to a quasi- periodic one at v — 0.087679 and there 
is transition to chaos below v = 0.087023. It is quite hard to discern the difference between a 
chaotic attractor and a quasi-periodic one by inspection of the signal of E(t). The return map is 
a. valuable tool in doing this, since the return map (of the energy minima for instance) of a quasi 
periodic attractor consists of lines with no foliations; an analogy is the return map obtained from a 
Poincare section of motion on a torus, in which case it is a circle whose perimeter becomes dense as 
the iterations increase. In the chaotic regime, the return map contains foldings and self-similarity 
as exhibited in several examples earlier. This type of chaotic motion is differentiated from what 
we term as unrecognizable chaos in that the tools employed by us to analyze the dynamics, yield 
information about the structure of the attractor. 

As can be seen from Figure 12, among the regions of different unsteady and complicated spatio- 
temporal motions there emerge stable fixed point attractors with increasingly higher modal be- 
havior, i.e. with shorter periods. In Figure 13 we give a summary of the computed attractors so 
far, as v decreases to approximately 0.006. The numbers on the figure denote the modal form of 
the solution, for example 3M is a trimodal attractor with the only non-zero Fourier coefficients 
being multiples of three. As already noted in the Appendix, one can construct steady states of 
increasingly shorter periods (i.e. higher modal behavior) by application of the similarity property 
(26). This construction does not imply stability but our numerical results produced by solving 
initial value problems, show that such fixed points are observable as large time solutions, given 
appropriate initial conditions. In what follows we use numerically computed solutions from each 
of the windows which support multimodal steady states (i.e. 3M, 4M, 5M, 6M, 7M, 8M and 9M) 
and show that they come from the unimodal fixed point attractor in interval A (see Figure 1). 

Making the assumption that the multimodal fixed points arise from the self-similar property 
(26), we first transform the ia windows which support these steady states by multiplying by N 2 
where N is the modal characterization of the attractor. This transformation gives the corresponding 
viscosity interval in the global fixed point attractor A of Figure 1. The results are summarized in 
Table 1. 
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Attractor 

v interval 

Transformed ^-interval (xJV 2 ) 

3M 

[.07465, .059922] 

[.67185,-539298] 

4M 

[.043848,-034625] 

[.701568, .554] 

5M 

[.026984,-022604] 

[.6746, .5651] 

6M 

[.01829, .015845] 

[.65844,.57042] 

7M 

[-012994,-01174] 

[.636706,-57526] 

8M 

[.009905,-009204] 

[.63392, .589056] 

9M 

[.007479, .007305] 

[.605799,-591705] 

10M 

.006 

.6 1 


Table 1 

Table 1 indicates that the transformed intervals for each of the multimodal attractors listed lie in 
the ^-interval which supports a unimodal global fixed point. In what follows we provide numerical 
evidence that each of the multimodal attractors at decreasing values of v, are indeed obtainable 
from the transformation (26). This is done for the profile corresponding to v = .6 in interval A, since 
this is seen to be common to all entries in the Table. For A r = 3, . . . , 10 we compute the steady state 
attractor at the corresponding values v = 0.6/iV 2 and denote the profile by u N (x). These profiles 
are shown in Figure 14; the period of corresponding profiles is seen to be 27r/A and the amplitude 
increases with N. Note that due to translation invariance each of the profiles can be shifted by a 
constant along the x-axis and it remains a periodic solution. This observation is particularly useful 
in the normalization construction that follows. Each profile is then normalized according to (26), 
that is we define a new function, U^(y) say, on [0,27 t] by applying the transformation 

Un(u) - j^u n (x/N). (25) 

This transformation was first noted by Frisch et al. [12] who went on to carry out the linear 
stability of these profiles for large N - see below. If the self-similarity property (26) is valid for 
these multimodal fixed points, then the Uiv(y) for N = 1 , . . ., 10 are identical to within a horizontal 
shift and equal to the profile obtained at u = 0.6. It was found that the computed profiles 4M-10M 
had to be shifted by 7 r in order to put them in phase with the globally attracting profile at v = 0.6 
(see earlier comments). The results are given in Figure 15 with different symbols used to denote 
different values of N. The solid line corresponds to the generating fixed point at v = 0.6. It is 
seen from the Figure that the self- similarity property is indeed the controlling mechanism in the 
generation of these high-modal steady states as v decreases. 

The data of Table 1 indicate the parameter range of the fixed point unimodal attractor which 
needs to be considered in order to transform, according to the similarity property, to the interval 
obtained numerically for a given multimodal attractor. In the stability studies of Frisch et al [12] 
and Papageorgiou et al [26], the underlying basic profiles are the multimodal ones indicated above 
and constructed in the Appendix, as the value of N or p increases. Such solutions are found to be 
stable in the i/-interval (0.6, 0.7), approximately. This result would imply that given a value of v 
in the stability window, the multimodal profiles obtained by the self similarity property (26), i.e. 
the one-parameter family Nu(Nx; v/N 2 ) for N = 1,2, ... as N gets large, are linearly stable. Our 
numerical results are in agreement with this theoretical prediction for N as large as 10. Note that 
computation of steady states with N > 10 has not been achieved yet because the window which 
supports such states has size of the order of machine precision. This may be achievable with higher 
precision arithmetic but such explorations are not undertaken here. 
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4 Conclusions 


We have given many quantitative features of solutions of the Kuramoto-Sivashinsky equations 
computed with general initial data and spatially periodic boundary conditions. Features such as 
strange attractors and periodic or chaotic bursting phenomena have been elucidated. Of particular 
interest is a set of computed multimodal solutions found at decreasing values of v and supported 
on smaller and smaller windows. It has been shown that these solutions (the last computed one is 
a decamodal profile, i.e. one with all Fourier coefficients zero except those which are multiples of 
10) derive from the unimodal fixed point attractor in .3 < v < 1, by a self- similar transformation 
property of the equation. It can be concluded that such solutions are at least linearly stable (in 
the sense that they are computable as large time solutions of a nonlinear initial value problem) 
and our numerical results are in full agreement with the stability theory of Frisch et al [12] and 
the modulation theory of Papageorgiou et al [26] who study the stability of such multimodal 
steady states. The numerical results given here provide strong support for the modulation stability 
theories. 


APPENDIX 

The following property of (2) is useful. It has been established numerically that for 0.307602 < 
v < 1 (this is interval A described in Figure 1 and Section 3) the KS has a global unimodal fixed 
point attractor. That is given fairly general initial conditions (for example a white noise spectrum) 
solutions evolve to a steady state with period 27T. Denote such steady states by U(x; v) with v in 
the interval A. For any constants p and c, equation (2) has the following two-parameter family of 
solutions 


u{x,t) = pU(p(x - ct);p 2 v) + c. (26) 

This follows by the scale and Galilean invariance of KS. For stationary waves c = 0, while peri- 
odicity implies that p is an integer. Property (26) can be used to generate steady states of (2) 
at geometrically decreasing values of the viscosity. For example if the 27r-periodic steady state is 
known at v = 0.8, we can construct a steady state at v = 0.8/2 2 by choosing p = 2: this solution 
has period 7r, i.e. is bimodal, and amplitude twice that of the starting steady state; this folding and 
scaling process can be repeated ad infinitum to obtain multimodal steady states at geometrically 
decreasing values of v. Clearly in the small v limit the amplitude of such steady states is of the 
order of v~ x ! 2 , which is also the order of the X 2 -norm of these steady state solutions. Even though 
the steady states exist and can be constructed as shown above, they are not necessarily stable and 
may not be observed in numerical experiments. What has been found, however, is that there is a 
subwindow of interval A which can be used to generate multimodal steady states which are linearly 
stable. 
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Figure 1 Schematic of the various attractors for larger values of v. The diagram is to scale 
and the various most strongly attracting solution regions are indicated. 
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Figure 3 Interval C: Periodic homoclinic bursts for v = 0.2. Time evolution of the 1st 
and 2nd Fourier modes. The graphs show the evolution of the coefficients of cos(ar), sin(z) and 
cos(2*), sin(2z). Higher mode amplitudes have much lower amplitudes. 
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Figure 4 Interval E: Period doubling route to chaos according to a Feigenbaum scenario. The 
Figure gives the phase planes of the energy as v decreases. The values of v and the number of 
relative minima in each curve are indicated on the Figure. Chaos sets in as the phase plane attains 
increasingly more turns which become infinitely many as an accumulation point is reached. 
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Figure 5 Interval E\: Feigenbaum chaos, v - 0.12118. Strange attractor produced by the 
return map of the energy minima. The Figure gives successive enlargements near the point ‘e’ 
marked on each subplot indicating the self-similar nature of the attractor. An overall enlargement 
of 3 x 10 6 takes place in going from 5(i) to 5(vii). 
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Figure 6 Interval F: Chaotic homoclinic bursts, v = 0.12115. The energy and its phase plane 
are shown. The bursts occur at regular intervals and are chaotic in nature. The steady attractor 
between bursts corresponds to the unimodal global steady attractor in interval A having v = 0.4846 
and the similarity transformation (26) (see text). 
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Figure 7 Interval F: Chaotic homoclinic bursts, v = 0.12115. Time history of the critical points 
of the energy (both maxima and minima). The absence of points indicates a constant energy. An 
extended chaotic region is seen between approximately t - 1000 - 1400, which then gets attracted 
to chaotic homoclinic bursts. The lower figure shows the return map of the energv minima. A 
strange attractor emerges and the point between bursts is indicated by an open circle. 
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Figure 8 Interval F: Chaotic homoclinic bursts, v = 0.12115. Self- similarity of the strange 
attractor. Successive enlargements are shown in the neighborhood of ‘o’. An overall enlargement 
of the order of 10 s takes place in going from 8(i) to 8(vi). 
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Figure 10 Interval F: Chaotic homoclinic bursts, v - 0.1. Return map of the energy minima 
indicating a strange attractor. The point P indicates the position in between bursts. 
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Figure 11 Interval F: Chaotic homoclinic bursts, v = 0.1. Evolution of the first 4 Fourier 
amplitudes together with the energy (sum of the Fourier amplitudes). It is seen that in between 
bursts the even Fourier modes are non-zero while the odd ones are zero, and the steady attractor 
is bimodal, then. During bursting an energy exchange takes place with odd modes gaining energy 
from even modes. The interval between bursts remains fairly constant over long time periods. 
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Figure 12 Schematic of the various computed attractors for smaller values of nu (u < 0.09). 
The diagram is not drawn to scale but window boundaries are given along with a brief description 
of the type of computed dynamics. The attractors are not global. 




nu 


Figure 13 Distribution of computed steady attractors for v < 0.09. The diagram, drawn to 
scale, shows the ^-intervals where steady state solutions were found. The notation 1M, 3M, etc. 
denotes profiles which are unimodal, trimodal etc.. The windows of support of higher modal steady 
states get decreasingly small as v decreases and their computation becomes prohibitively difficult. 
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Figure 14 Computed profiles corresponding to values of u taken from the global attractor 
for large v (interval A in Figure 1) and the intervals 3M, 4M, 5M, 6M, 7M, 8M, 9M and 10M. 
The corresponding values of nu are indicated on the Figures and are chosen in order to check the 
self-similarity of the solutions (see text). 
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Figure 15 Similarity of steady states as v decreases. Collapse of the profiles given in Figure 
14 onto the universal profile corresponding to u = 0.6. Different symbols correspond to data taken 
from the indicated intervals. The similarity property (26) is used in the construction of the Figure 
- see text for details. 
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